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Abstract. 

In  this  report  we  list  the  Fortran  subroutines  for 
evaluating  the  confluent  hypergeomet r i c  functions  M(a,b;x)  and 
U(a,b;x) .   These  subroutines  use  the  stable  recurrence 
relations  given  e.g.  in  Wimp. 
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Int roduct ion 

It  is  well  known  that  the  ordinary  differential  equation 


d  y     ,  *.        n  d  y  „ 

X  ^2  +  C1-)  d^  "  ay  =  ° 


has  a  solution 

y(x)  =  AM(a,l;x)  +  BU(a,l;x) 

if  a  is  not  a  negative  integer. 

This  problem  arises  e.g.  when  solving  the  linearized 
shallow  water  equations  with  the  full  linear  variation  in  depth 
included  (see  Williams,  Staniforth  and  Neta,   [1] ) . 

The  computation  of  the  confluent  hypergeomet r i c  functions 
i  s  based  on  the  Miller  a  Igor  ithm  (see  e.g.  Wimp,   [2]).    In 
general,  one  has  a  second  order  difference  equation 

z(n)  +  a(n)z(n+l)  +  b(n)z(n+2)  =0,   n  >  0,   b(n)  ^  0  . 

If  b(n)  =  0  for  some  n,  in  some  cases  one  can  make  a  change  of 
variable  Y(n)  =  A(n)z(n)  which  will  produce  an  equation  of  the 
desired  type.   Let  w(n)  be  a  nontrivial  solution  and  the  sum  of 

the  n  /rmalizing  series 


S  =  §   c(k)w(k)  #  0 
k=0 


is  known.    Let  N  bo  a   large  i nt c^cr    and  define  zN(n),  0  <  n  < 

N+l  .  by 


zN  (  n  ) 


0    n  =  N+l 


1    n  =  N 


zN(n)  +  a(n)zN(n  +  l)  +  b(n)zN(n  +  2)  =0,     n  =  N-l 1,  0 


One  can  approximate  w(n)  by  wN(n) 


(n)  =  SzN(n)/S 


wnere 


N 
SN  ~  T,     c(k)zN(k) 
k=0 


The  algorithm  is  said  to  converge  if 


vN  (n )   —   w  ( n )   as   N  — »  oo 


The  function  M(a,b:x)  satisfies  the  recurrence  relat ion 


(2n+b+2)(n+a)z(n)     (2n+b+l ) { ( 2a-b)  + 


(2n+b) (2n+b+2) 


}z(n+l) 


(2n+b) (n+b+l-a)z(n+2)  =  0  . 


Tiie    minimal     solution 


w(n)     =    "      \    >u    M(a+n,2n+b;x) 

(  b  )  2n 


where 


r(n+c) 

Cc;n  -     r(c) 


The  normalization  relationship  used  in  our  subroutine  is 


S  =  b-1  =  §   (  */  (b-l)k(b+2k-l)w(k) 

k=0  k- 


An  obvious  modification  must  be  made  if  b  =  1.   The  algorithm 
is  not  defined  if  b.  b+l-a,  a  are  negative  integers  or  zero. 
The  function  U(a,b;x)  satisfies  the  relationship 

(n+a) (n+a+l-b)z(n)  -  ( n  +  1 )  [2 (n+a+1 ) +x-b] z (n  +  1 ) 

+  (n+1) (n+2)z(n+2)  =  0  . 

The  minimal  solution  is 


.  N     xn(a)n(a+l-b)n 
w  (  n  )  = — ^-^ —    U  (  a+n  .  b  :  x  ) 


for  |  arg  x|  <  tt  .        A  normalization  relation  is 


oo 
1  =  £   w(k) 
k=0 


In  the  next  section  we  give  a  listing  of  the  Fortran 
subrout  i  nes . 


Subroutine  Miller 

SUBROUTINE  MILLER(N , ALPHA , BETA . X , S . SS . COEFF) 
INTEGER  N 

REAL*8   ALPHA, BETA, X , SS 
REAL*8   S (0:1000) 
EXTERNAL  COEFF 
C     USES  THE  J.C.P.  MILLER  ALGORITHM  TO  COMPUTE 
C     S(0:N). 
C     BEGIN 

INTEGER  NN.K 
REAL*8  T.D.EPS,A,B,C 
REAL*8  0LDS(0:1000) 
EPS  =  0.000000001 
C         INITIALIZE  OLDS. 
DO  1  K  =  0,  1000 
OLDS(K)  =  0 

1  CONTINUE 

C        CHOOSE  INITIAL  NN . 

NN  =  N  +  2 
C         INITIALIZE  K,  S  AND  T. 

2  K  =  NN 
S(K+1 )  =  0 
S(K)    =1 

CALL  COEFF (K , ALPHA . BETA , X , A , B , C) 
T       =  2*C*S(K) 
C         TAKE  A  BACKWARD  RECURRENCE  STEP  AND  UPDATE  IT 

3  K  =  K  -  1 

CALL  COEFF (K . ALPHA , BETA , X . A . B , C) 
S(K)  =  A*S(K+1)  +  B*S(K+2) 
C         CHECK  FOR  OVERFLOW  AND  RESCALE  IF  NECESSARY. 
D=  DABS(S(K) ) 
IF  (D  .GT.  1.D30)  THEN 
C         BEGIN 

CALL  SCALE (K,NN,S.T.D) 
END  IF 

IF  (K  .GT.  0)  THEN 
C        BEG  I N 

T  =  T  +  2*C*S(K) 

GO  TO  3 
END  IF 

T  =  T  +  C*S(0) 
DO  4  K  =  0 ,  N 

S(K)  =  S(K)/T 

4  CONTINUE 

C         TEMPORARY  PRINT  STATEMENT. 

C        PRINT*,  S(0) 

C        TEST  FOR  CONVERGENCE. 

D  =  0 

DO  5  K  =  0,  N 

5  D  =  D  +  S(K)**2 
CONTINUE 

D  =  DSQRT(D) 
T  =  0 


DO  6  K  =  0,  N 

T  =  T  +  (S(K)  -  0LDS(K))**2 

6  CONTINUE 

T  =  DSQRT(T) 
C        TAKE  ANOTHER  STEP  IF  NO  CONVERGENCE. 

IF  (T  .GT.  EPS*D)  THEN 
C         BEGIN 

NN  =  2*NN 

DO  7  K  =  0,  N 

OLDS(K)  =  S(K) 

7  CONTINUE 

IF(NN  .LE.  1000)  GO  TO  2 
PRINT  999, NN, ALPHA, BETA, X,T 
999      FORMAT  ('  **  NO  CONVERGENCE  **  NN  AP  CP  X  T  ',I5,4E14.7) 
END  IF 
SS=S(0) 
RETURN 
END 


SUBROUT I NE  COEFF ( N , ALPHA , BETA , X , A , B . C) 

INTEGER  N 

REAL*8  ALPHA. BETA, X, A, B,C 
C     COMPUTES  COEFFICIENTS  USED  BY  .)  .  C .  P  .  MILLER  ALGORITHM  FOR 
C     A  CONFLUENT  HYPERGEOMETRIC  FUNCTION   M(a,b:x) 
C     SEE  JET  WIMP.  COMPUTATION  WITH  RECURRENCE  RELATIONS, 
C     PITMAN  1984  PP.  61-62 
C     BEGIN 

INTEGER  M.K 
REAL*g.T,U,V,W 


C 


S 

=  2*ALPHA  - 

BETA 

T 

=  N  +  ALPHA 

M 

=  2*N 

U 

=  M  +  BETA 

V 

=  U  +  1 

w 

=  V  +  1 

A 

=  (S/W  +  U/X)*V/T 

B 

=  (N  +  BETA 

-  ALPHA 

T 

=  1 

IF 

1  (N  .GT.  0) 

THEN 

BEG  I N 

S  =  BETA  - 

1 

DO  1  K  =  1 

,  N-l 

T  =  -T* 

(1+S/K) 

CONTINUE 

T  =  -T*(l+S/M) 

END 

C 

=  T 

RETURN 

+  1)*U/T/W 


SUBROUT I NE  SCALE ( K . N , S . T , D ) 
INTEGER  N.K 
REALMS  T.D 
REAL*8  S (0:1000) 
C     BEG  I N 

INTEGER  J 

T  =  T/D 

DO  1  J  =  K,  N 

S(J)  =  S(J)/D 
1     CONTINUE 
RETURN 
END 


SUBROUTI NE  COEFU (N , ALPHA , BETA , X , A , B , C) 
INTEGER  N 

REAL*8  ALPHA,  BETA , X , A , B , C 
C     COMPUTES  COEFFICIENTS  USED  BY  J.C.P.  MILLER  ALGORITHM  FOR 
C     A  CONFLUENT  HYPERGEOMETRIC  FUNCTION   U(a,b;x) 
C     SEE  JET  WIMP,  COMPUTATION  WITH  RECURRENCE  RELATIONS, 
C     PITMAN  1984  PP.  63-64 
C     BEGIN 

INTEGER  M.K 

REAL*8   S,T,U,V,W 

S  =    ALPHA  +  QFLOAT(N) 

T  =  S  +  1 .DO 

U  =  S*(T  -  BETA) 

V  =  QFLOAT(N  +  1) 

W  =  V  +  1 . DO 

A  =  (2*T  +  X  -  BETA)*V/U 

B  =  -  V*W/U 

C  =  1 

RETURN 
END 


Remark:   The  program  "that:  calls  Miller  must:  supply  as  a  last 
parameter  either  COEFF  (for  M)  or  COEFU  (for  U) . 


The  subroutines  are  available  on  a  diskette  from  either  author 
upon  request.   These  subroutines  were  tested  extensively  for 
various  values  of  a,  b  and  x. 

Remark :    If  the  parameter  is  a  negative  integer,  the  solul  ion 
of  the  differential  equation  is 


y  =  ALn(x)  +  B{ln|x|Ln(x)  +  g  /3mxm} 

m=0 


where  n  =  -a. 

Ln(x)  are  Laguerre  polynomials  whose  coefficients  fi| 

s  at  i  s  f  y 

_  i  -  n  - 1  :    o         n 

l 

q  j  =  -  n 

The    coefficients    i3m    satisfy 

/            ^   d              ( ->           2  (  m  -  n  )           \ 
(m-n)iJm    +    V     ~         m+1         °m 
/^rn  +  1     =     7~VT2 m    =     1 n_1 


(m+1) 


0m    =     ~ 5    «n  m    =    n 

(n+1)2 


0m         D1^1    /?m.i  ™    =    n  +  1  ,  n+2 

m 
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